1//////////////////////////////////////////////////////////////////////
2// LibFile: convex_hull.scad
3// Functions to create 2D and 3D convex hulls.
4// To use, add the following line to the beginning of your file:
5// ```
6// include <BOSL/convex_hull.scad>
7// ```
8// Derived from Linde's Hull:
9// - https://github.com/openscad/scad-utils
10//////////////////////////////////////////////////////////////////////
11
12include <BOSL/math.scad>
13
14
15
16// Section: Generalized Hull
17
18// Function: convex_hull()
19// Usage:
20// convex_hull(points)
21// Description:
22// When given a list of 3D points, returns a list of faces for
23// the minimal convex hull polyhedron of those points. Each face
24// is a list of indexes into `points`.
25// When given a list of 2D points, or 3D points that are all
26// coplanar, returns a list of indices into `points` for the path
27// that forms the minimal convex hull polygon of those points.
28// Arguments:
29// points = The list of points to find the minimal convex hull of.
30function convex_hull(points) =
31 !(len(points) > 0) ? [] :
32 len(points[0]) == 2 ? convex_hull2d(points) :
33 len(points[0]) == 3 ? convex_hull3d(points) : [];
34
35
36
37// Section: 2D Hull
38
39// Function: convex_hull2d()
40// Usage:
41// convex_hull2d(points)
42// Description:
43// Takes a list of arbitrary 2D points, and finds the minimal convex
44// hull polygon to enclose them. Returns a path as a list of indices
45// into `points`.
46function convex_hull2d(points) =
47 (len(points) < 3)? [] : let(
48 a=0, b=1,
49 c = _find_first_noncollinear([a,b], points, 2)
50 ) (c == len(points))? _convex_hull_collinear(points) : let(
51 remaining = [ for (i = [2:len(points)-1]) if (i != c) i ],
52 ccw = triangle_area2d(points[a], points[b], points[c]) > 0,
53 polygon = ccw? [a,b,c] : [a,c,b]
54 ) _convex_hull_iterative_2d(points, polygon, remaining);
55
56
57// Adds the remaining points one by one to the convex hull
58function _convex_hull_iterative_2d(points, polygon, remaining, _i=0) =
59 (_i >= len(remaining))? polygon : let (
60 // pick a point
61 i = remaining[_i],
62 // find the segments that are in conflict with the point (point not inside)
63 conflicts = _find_conflicting_segments(points, polygon, points[i])
64 // no conflicts, skip point and move on
65 ) (len(conflicts) == 0)? _convex_hull_iterative_2d(points, polygon, remaining, _i+1) : let(
66 // find the first conflicting segment and the first not conflicting
67 // conflict will be sorted, if not wrapping around, do it the easy way
68 polygon = _remove_conflicts_and_insert_point(polygon, conflicts, i)
69 ) _convex_hull_iterative_2d(points, polygon, remaining, _i+1);
70
71
72function _find_first_noncollinear(line, points, i) =
73 (i>=len(points) || !collinear_indexed(points, line[0], line[1], i))? i :
74 _find_first_noncollinear(line, points, i+1);
75
76
77function _find_conflicting_segments(points, polygon, point) = [
78 for (i = [0:len(polygon)-1]) let(
79 j = (i+1) % len(polygon),
80 p1 = points[polygon[i]],
81 p2 = points[polygon[j]],
82 area = triangle_area2d(p1, p2, point)
83 ) if (area < 0) i
84];
85
86
87// remove the conflicting segments from the polygon
88function _remove_conflicts_and_insert_point(polygon, conflicts, point) =
89 (conflicts[0] == 0)? let(
90 nonconflicting = [ for(i = [0:len(polygon)-1]) if (!in_list(i, conflicts)) i ],
91 new_indices = concat(nonconflicting, (nonconflicting[len(nonconflicting)-1]+1) % len(polygon)),
92 polygon = concat([ for (i = new_indices) polygon[i] ], point)
93 ) polygon : let(
94 before_conflicts = [ for(i = [0:min(conflicts)]) polygon[i] ],
95 after_conflicts = (max(conflicts) >= (len(polygon)-1))? [] : [ for(i = [max(conflicts)+1:len(polygon)-1]) polygon[i] ],
96 polygon = concat(before_conflicts, point, after_conflicts)
97 ) polygon;
98
99
100
101// Section: 3D Hull
102
103// Function: convex_hull3d()
104// Usage:
105// convex_hull3d(points)
106// Description:
107// Takes a list of arbitrary 3D points, and finds the minimal convex
108// hull polyhedron to enclose them. Returns a list of faces, where
109// each face is a list of indexes into the given `points` list.
110// If all points passed to it are coplanar, then the return is the
111// list of indices of points forming the minimal convex hull polygon.
112function convex_hull3d(points) =
113 (len(points) < 3)? list_range(len(points)) : let (
114 // start with a single triangle
115 a=0, b=1, c=2,
116 plane = plane3pt_indexed(points, a, b, c),
117 d = _find_first_noncoplanar(plane, points, 3)
118 ) (d == len(points))? /* all coplanar*/ let (
119 pts2d = [ for (p = points) xyz_to_planar(p, points[a], points[b], points[c]) ],
120 hull2d = convex_hull2d(pts2d)
121 ) hull2d : let(
122 remaining = [for (i = [3:len(points)-1]) if (i != d) i],
123 // Build an initial tetrahedron.
124 // Swap b, c if d is in front of triangle t.
125 ifop = in_front_of_plane(plane, points[d]),
126 bc = ifop? [c,b] : [b,c],
127 b = bc[0],
128 c = bc[1],
129 triangles = [
130 [a,b,c],
131 [d,b,a],
132 [c,d,a],
133 [b,d,c]
134 ],
135 // calculate the plane equations
136 planes = [ for (t = triangles) plane3pt_indexed(points, t[0], t[1], t[2]) ]
137 ) _convex_hull_iterative(points, triangles, planes, remaining);
138
139
140// Adds the remaining points one by one to the convex hull
141function _convex_hull_iterative(points, triangles, planes, remaining, _i=0) =
142 _i >= len(remaining) ? triangles :
143 let (
144 // pick a point
145 i = remaining[_i],
146 // find the triangles that are in conflict with the point (point not inside)
147 conflicts = _find_conflicts(points[i], planes),
148 // for all triangles that are in conflict, collect their halfedges
149 halfedges = [
150 for(c = conflicts, i = [0:2]) let(
151 j = (i+1)%3
152 ) [triangles[c][i], triangles[c][j]]
153 ],
154 // find the outer perimeter of the set of conflicting triangles
155 horizon = _remove_internal_edges(halfedges),
156 // generate a new triangle for each horizon halfedge together with the picked point i
157 new_triangles = [ for (h = horizon) concat(h,i) ],
158 // calculate the corresponding plane equations
159 new_planes = [ for (t = new_triangles) plane3pt_indexed(points, t[0], t[1], t[2]) ]
160 ) _convex_hull_iterative(
161 points,
162 // remove the conflicting triangles and add the new ones
163 concat(list_remove(triangles, conflicts), new_triangles),
164 concat(list_remove(planes, conflicts), new_planes),
165 remaining,
166 _i+1
167 );
168
169
170function _convex_hull_collinear(points) =
171 let(
172 a = points[0],
173 n = points[1] - a,
174 points1d = [ for(p = points) (p-a)*n ],
175 min_i = min_index(points1d),
176 max_i = max_index(points1d)
177 ) [min_i, max_i];
178
179
180
181function _remove_internal_edges(halfedges) = [
182 for (h = halfedges)
183 if (!in_list(reverse(h), halfedges))
184 h
185];
186
187
188function _find_conflicts(point, planes) = [
189 for (i = [0:len(planes)-1])
190 if (in_front_of_plane(planes[i], point))
191 i
192];
193
194
195function _find_first_noncoplanar(plane, points, i) =
196 (i >= len(points) || !coplanar(plane, points[i]))? i :
197 _find_first_noncoplanar(plane, points, i+1);
198
199
200// vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap